library(tidyverse)
library(openxlsx) #read excel file
library(ggrepel) #label text
library(vegan) #Biodiversity Index
library(sf) #Map
library(scales)
library(RColorBrewer)
library(treemapify)
source("Custom_function.R")
Read Nui Chua data with read.xlsx() function
#library(openxlsx)
filename <- "Plot_NC-BD_20230925_v4.xlsx"
# NuiChua plot
NCplot <- read.xlsx(filename, sheet = "NC_census2")
#Bidoup - NuiBa plot
#BDplot <- read.xlsx(filename, sheet = "BDNB_census2")
Display the Structure of data and view some first lines.
dbh: Diametre in mm;
lx,ly: tree coordinates to zero position of quadrat (20x20m);
gx,gy: tree coordinates to zero position of plot (100x100m).
str(NCplot)
## 'data.frame': 17383 obs. of 17 variables:
## $ PlotID : num 1 1 1 1 1 1 1 1 1 1 ...
## $ census : chr "C2" "C2" "C2" "C2" ...
## $ tag : chr "01-0048" "01-0049" "01-0049" "01-0050" ...
## $ stemtag : chr "01-0048_0" "01-0049_0" "01-0049_1" "01-0050_0" ...
## $ spcode : chr "dimosp2" "crotth" "crotth" "dryppo" ...
## $ quadrat : chr "0000" "0000" "0000" "0000" ...
## $ lx : num 1.119 1.545 1.545 1.188 0.327 ...
## $ ly : num 1.54 3.21 3.21 3.44 3.73 ...
## $ gx : num 1.119 1.545 1.545 1.188 0.327 ...
## $ gy : num 1.54 3.21 3.21 3.44 3.73 ...
## $ dbh : num 61 41 28 69 160 ...
## $ codes : chr NA "TNBT" "TNBT" "TBV" ...
## $ hom : chr NA NA NA NA ...
## $ date : chr "2021-03-18" "2021-03-18" "2021-03-18" "2021-03-18" ...
## $ family : chr "Euphorbiaceae" "Euphorbiaceae" "Euphorbiaceae" "Putranjivaceae" ...
## $ genus : chr "Dimorphocalyx" "Croton" "Croton" "Drypetes" ...
## $ scientificName: chr "Dimorphocalyx sp.2" "Croton thorelii" "Croton thorelii" "Drypetes poilanei" ...
head(NCplot)
## PlotID census tag stemtag spcode quadrat lx ly gx
## 1 1 C2 01-0048 01-0048_0 dimosp2 0000 1.118787 1.544663 1.118787
## 2 1 C2 01-0049 01-0049_0 crotth 0000 1.544927 3.205741 1.544927
## 3 1 C2 01-0049 01-0049_1 crotth 0000 1.544927 3.205741 1.544927
## 4 1 C2 01-0050 01-0050_0 dryppo 0000 1.188361 3.440553 1.188361
## 5 1 C2 01-0051 01-0051_0 syzysp2 0000 0.327383 3.727546 0.327383
## 6 1 C2 01-0052 01-0052_0 dimosp2 0000 1.127484 4.179776 1.127484
## gy dbh codes hom date family genus
## 1 1.544663 61.0 <NA> <NA> 2021-03-18 Euphorbiaceae Dimorphocalyx
## 2 3.205741 41.0 TNBT <NA> 2021-03-18 Euphorbiaceae Croton
## 3 3.205741 28.0 TNBT <NA> 2021-03-18 Euphorbiaceae Croton
## 4 3.440553 69.0 TBV <NA> 2021-03-18 Putranjivaceae Drypetes
## 5 3.727546 160.5 <NA> <NA> 2021-03-18 Myrtaceae Syzygium
## 6 4.179776 148.5 <NA> <NA> 2021-03-18 Euphorbiaceae Dimorphocalyx
## scientificName
## 1 Dimorphocalyx sp.2
## 2 Croton thorelii
## 3 Croton thorelii
## 4 Drypetes poilanei
## 5 Syzygium sp.2
## 6 Dimorphocalyx sp.2
Use n_distinct() counts the number of distinct value in
column: family, genus, scientificName
NCplot |>
summarise("family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName))
## family genus species
## 1 45 76 93
(spByFam <- NCplot |>
group_by(family) |>
summarise(n = n_distinct(spcode)) |>
arrange(desc(n)))
## # A tibble: 45 Ă 2
## family n
## <chr> <int>
## 1 Euphorbiaceae 7
## 2 Rubiaceae 7
## 3 Phyllanthaceae 6
## 4 Annonaceae 5
## 5 Ebenaceae 5
## 6 Fabaceae 4
## 7 Lamiaceae 4
## 8 Celastraceae 3
## 9 Melastomataceae 3
## 10 Myrtaceae 3
## # âč 35 more rows
Figure in treemap layout
library(treemapify)
ggplot(spByFam, aes(
area = n,
fill = n,
label = paste(family, n, sep = "\n")
)) +
geom_treemap(start = "bottomleft") +
geom_treemap_text(place = "centre") +
scale_fill_distiller(palette = "Greens", direction = 1) +
labs(fill = "Species")
Species accumulation curves is a graph recording the cumulative number of species recorded in a particular environment. It used to estimate the number of species in a particular area.
Here, we use b_specaccum() function.
curveNC <- b_specaccum(NCplot)
head(curveNC)
## subplot richness S PlotName
## 1 1 35.24000 400
## 2 2 48.11000 800
## 3 3 56.08174 1200
## 4 4 61.81360 1600
## 5 5 66.23241 2000
## 6 6 69.78524 2400
#FIGURE
ggplot(curveNC, aes(x = S, y = richness)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
geom_text(
aes(label = round(richness, 1)),
size = 4,
vjust = -0.5,
check_overlap = T
) +
theme_bw(16) +
labs(x = "Survey area (m2)", y = "Number of species")
We usually group trees into 3 size-class of DBH:
small tree [1-5)cm;
medium tree [5-20)cm;
big tree >=20 cm
NCplot <- NCplot |>
mutate(
dbh_cm = dbh / 10,
# dbh value in mm
gDBH = case_when(dbh_cm >= 20 ~ ">=20",
dbh_cm >= 5 ~ "[5-20)",
TRUE ~ "[1-5)")
)
head(NCplot)
We count the number of stems at each diameter class.
(stemByDBH <- NCplot |>
count(gDBH) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n)))
## gDBH n freq
## 1 [1-5) 14221 0.818
## 2 [5-20) 2951 0.170
## 3 >=20 211 0.012
Figure
stemByDBH |>
ggplot(aes(x = gDBH, y = n)) +
geom_col(aes(fill = n)) +
geom_text(aes(label = paste0(n, " (", percent(freq), ")")), vjust = -0.5) +
scale_fill_viridis_c() +
theme_bw(14) +
labs(x = "DBH (cm)", y = "Nb of stem", fill = "Nb of stem")
We can also count the number of species at each diameter class.
(
SpByDBH <- NCplot |>
group_by(gDBH) |>
summarise(
"family" = n_distinct(family),
"genus" = n_distinct(genus),
"species" = n_distinct(scientificName)
)
)
## # A tibble: 3 Ă 4
## gDBH family genus species
## <chr> <int> <int> <int>
## 1 >=20 26 33 39
## 2 [1-5) 41 71 87
## 3 [5-20) 38 59 72
Figure
SpByDBH |>
tidyr::pivot_longer(-gDBH, names_to = "rank", values_to = "n") |>
ggplot(aes(x = gDBH, y = n, fill = rank)) +
geom_col(position = position_dodge()) +
geom_text(aes(x = gDBH, y = n, label = n),
position = position_dodge(width = 0.9),
vjust = -0.5) +
scale_fill_viridis_d() +
theme_bw(15) +
labs(x = "DBH class", y = "Number", fill = "Rank")
Draw a blank plot with draw_blank_plot_1ha()
function.
p <- draw_blank_plot_1ha()
p
Add add tree to map
p + geom_point(data = NCplot,
aes(
x = gx,
y = gy,
size = dbh_cm,
color = gDBH
),
alpha = 0.7) +
labs(color = "DBH class", size = "DBH(cm)")
For big tree
NCplotA <- NCplot |> filter(gDBH == ">=20")
p + geom_point(
data = NCplotA,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7,
shape = 21,
fill = "blue",
color = "white"
) +
labs(size = "DBH(cm)")
For medium tree
NCplotB <- NCplot |> filter(gDBH == "[5-20)")
p + geom_point(
data = NCplotB,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7,
shape = 21,
fill = "blue",
color = "white"
) +
labs(size = "DBH(cm)")
For small tree
NCplotC <- NCplot |> filter(gDBH == "[1-5)")
p + geom_point(
data = NCplotC,
aes(x = gx, y = gy, size = dbh_cm),
alpha = 0.7,
shape = 21,
fill = "blue",
color = "white"
) +
labs(size = "DBH(cm)")
For a particular species, px: Dialium cochinchinense
## So do phan bo toan bo cay trong o mau
NCplotSP <- NCplot |> filter(scientificName == "Dialium cochinchinense")
p + geom_point(data = NCplotSP,
aes(
x = gx,
y = gy,
size = dbh_cm,
color = gDBH
),
alpha = 0.7) +
labs(color = "DBH class", size = "DBH(cm)")
We can also try with a detail DBH class
NCplot <- NCplot |>
mutate(
dbh_cm = dbh / 10 , # dbh in mm
gDBH2 = case_when(
dbh_cm >= 100 ~ ">=100",
dbh_cm >= 90 ~ "[90-100)",
dbh_cm >= 80 ~ "[80-90)",
dbh_cm >= 70 ~ "[70-80)",
dbh_cm >= 60 ~ "[60-70)",
dbh_cm >= 50 ~ "[50-60)",
dbh_cm >= 40 ~ "[40-50)",
dbh_cm >= 30 ~ "[30-40)",
dbh_cm >= 20 ~ "[20-30)",
dbh_cm >= 10 ~ "[10-20)",
TRUE ~ "[1-10)")
)
stemByDBH2 <- NCplot |>
count(gDBH2) |>
mutate(freq = round(n / sum(n), 3)) |>
arrange(desc(n))
stemByDBH2
## gDBH2 n freq
## 1 [1-10) 16210 0.933
## 2 [10-20) 962 0.055
## 3 [20-30) 162 0.009
## 4 [30-40) 33 0.002
## 5 [40-50) 6 0.000
## 6 [50-60) 6 0.000
## 7 [60-70) 2 0.000
## 8 >=100 1 0.000
## 9 [70-80) 1 0.000
stemByDBH2 |>
ggplot(aes(x = gDBH2, y = n)) +
geom_col() +
geom_text(aes(label = n), vjust = -0.5) +
#scale_fill_manual(values = c("grey","white"))+
theme_bw(14) +
labs(x = "DBHClass (cm)", y = "Number of stem")
The Importance Value Index (IVI) in Ecology is the quantitative measure of how dominant a species is in a given ecosystem. IVI is calculated from 3 values:
relative density of each species (RD)
relative basal area of each species (RBA)
relative frequency of each species (this can be omitted in 1 plot)
\[IVI={\left(relative.density + relative.basal.area \right)}\]
Basal area (BA) of each species can be calculated from diameter measured: \(BA = pi * R^2\)
Relative basal area (RBA): \(RBA = BA / sum(BA)\)
Relative density (RD) = number of individuals of each species / total
NCplot <- NCplot |>
mutate(BA = ((dbh/2)/1000)^2*pi) # unit: m2)
IVI <- NCplot |>
group_by(scientificName) |>
summarise(n = n(), sBA = sum(BA)) |>
mutate(
RBA = sBA / sum(sBA) * 100,
RD = n / sum(n) * 100,
IVI = RBA + RD,
rank = rank(-IVI)
) |>
arrange(desc(IVI))
Figure
ggplot(IVI, aes(x = rank, y = IVI)) +
geom_line() +
geom_point() +
geom_label_repel(
data = IVI[1:5, ],
aes(label = paste0("italic('", scientificName, "')")),
fill = NA,
color = "black",
size = 5,
label.size = NA,
min.segment.length = 0,
segment.linetype = 2,
segment.size = 0.4,
parse = TRUE
) +
theme_bw(14) +
labs(x = "Species rank", y = "IVI (%)")
A diversity index is a quantitative measure that reflects how many different types (such as species) there are in a dataset (a community). These indices are statistical representations of biodiversity in different aspects (richness, evenness, and dominance).
Richness is the total number of unique species
As species richness and evenness increase, so diversity increases. The Shannon Diversity Index (H) is a way to measure the diversity of species in a community. The higher the value of H, the higher the diversity of species in a particular community. The lower the value of H, the lower the diversity. A value of H = 0 indicates a community that only has one species.
The Shannon Equitability Index is a way to measure the evenness of species in a community. The term âevennessâ simply refers to how similar the abundances of different species are in the community. This value ranges from 0 to 1 where 1 indicates complete evenness.
\[ Evenness = Shannon Index (H) / ln(richness) \]
Simpsonâs Index (D) measures the probability that two individuals randomly selected from a sample will belong to the same species
With this index, 0 represents infinite diversity and 1, no diversity. That is, the bigger the value of D, the lower the diversity. This is neither intuitive nor logical, so to get over this problem, D is often subtracted from 1 to give:
Simpsonâs Index of Diversity 1 - D
The value of this index also ranges between 0 and 1, but now, the greater the value, the greater the sample diversity. This makes more sense. In this case, the index represents the probability that two individuals randomly selected from a sample will belong to different species.
library(vegan)
# Simpson's Index of Diversity
D <- diversity(IVI$n, index = "simpson")
# Shannon
H <- diversity(IVI$n, index = "shannon")
# Richness
S <- specnumber(IVI$n)
# Evenness
E <- H/log(S)
# Create data frames
res <- data.frame("Richness"=S,"Simpson"=D,"Shannon"=H,"Evenness"=E)
res
## Richness Simpson Shannon Evenness
## 1 93 0.9014608 2.927733 0.645928
Create calcIVI() to calculator Importance Value Index
and figureIVI() to for plotting figure
calcIVI <- function(dataplot) {
IVI <- dataplot |>
mutate(G = ((dbh / 2) / 1000)^2 * pi) |>
group_by(scientificName) |>
summarise(n = n(), sG = sum(G)) |>
mutate(
UT = sG / sum(sG) * 100,
MD = n / sum(n) * 100,
IVI = UT + MD,
rank = rank(-IVI)
) |>
arrange(desc(IVI))
return(IVI)
}
figureIVI <- function(IVI) {
ggplot(IVI, aes(x = rank, y = IVI)) +
geom_line() +
geom_point() +
geom_label_repel(
data = IVI[1:5, ],
aes(label = paste0("italic('", scientificName, "')")),
fill = NA,
color = "black",
size = 5,
label.size = NA,
min.segment.length = 0,
segment.linetype = 2,
segment.size = 0.4,
parse = TRUE
) +
theme_bw(14) +
labs(x = "Species rank", y = "IVI (%)")
}
Create biodivIndex() to calculator Diversity Index
biodivIndex <- function(IVI) {
H <- diversity(IVI$n, index = "shannon")
D <- diversity(IVI$n, index = "simpson")
S <- specnumber(IVI$n)
E <- H / log(S)
res <- data.frame(
"Richness" = S,
"Simpson" = D,
"Shannon" = H,
"Evenness" = E
)
return(res)
}
Big tree (DBH >= 20 cm)
IVIA <- calcIVI(NCplotA)
head(IVIA)
## # A tibble: 6 Ă 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.1 24 3.20 22.8 11.4 34.1 1
## 2 Syzygium sp.2 37 2.05 14.6 17.5 32.1 2
## 3 Drypetes poilanei 24 1.05 7.48 11.4 18.9 3
## 4 Diospyros sp.2 21 1.19 8.44 9.95 18.4 4
## 5 Dialium cochinchinense 18 1.11 7.87 8.53 16.4 5
## 6 Dehaasia sp. 10 0.753 5.36 4.74 10.1 6
figureIVI(IVIA)
biodivIndex(IVIA)
## Richness Simpson Shannon Evenness
## 1 39 0.9170055 2.919805 0.7969853
Medium tree (DBH in 5-20 cm)
IVIB <- calcIVI(NCplotB)
head(IVIB)
## # A tibble: 6 Ă 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Syzygium sp.2 549 5.64 26.0 18.6 44.6 1
## 2 Drypetes poilanei 521 3.76 17.3 17.7 35.0 2
## 3 Diospyros sp.2 227 1.83 8.43 7.69 16.1 3
## 4 Cleistanthus sp. 327 1.08 4.97 11.1 16.1 4
## 5 Walsura bonii 136 0.868 4.00 4.61 8.61 5
## 6 Diospyros sp.1 116 0.853 3.93 3.93 7.86 6
figureIVI(IVIB)
biodivIndex(IVIB)
## Richness Simpson Shannon Evenness
## 1 72 0.9062302 2.945215 0.6886708
Small tree (DBH in 1-5 cm)
IVIC <- calcIVI(NCplotC)
head(IVIC)
## # A tibble: 6 Ă 7
## scientificName n sG UT MD IVI rank
## <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Cleistanthus sp. 2767 1.52 22.8 19.5 42.2 1
## 2 Dimorphocalyx sp.2 3215 1.10 16.5 22.6 39.1 2
## 3 Goniothalamus sp. 2230 1.13 17.0 15.7 32.6 3
## 4 Drypetes poilanei 599 0.465 6.97 4.21 11.2 4
## 5 Cleistanthus monoicus 555 0.340 5.09 3.90 8.99 5
## 6 Mallotus sp. 726 0.176 2.63 5.11 7.74 6
figureIVI(IVIC)
biodivIndex(IVIC)
## Richness Simpson Shannon Evenness
## 1 87 0.8767146 2.706991 0.6061457
Above-ground biomass (AGB)
Below-ground biomass (BGB)
Chave et al (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1): 87â99.
AGB=Total aboveground biomass (kg)Â
D=DBH (cm)Â
H=Height (m)Â
Ï=Wood density (g/cm3)Â
BA=Basal area (cm2)Â
| Equation | |
| Tropical forest DRY | \(AGB = exp(-2.187 + 0.916 * ln(Ï * D^2 * H)\) |
| \(AGB = Ï * exp(-0.667 + 1.784 * ln(D) + 0.207 * (ln(D))^2 â 0.0281 * (ln(D))^3)\) | |
| Tropical forest MOIST | \(AGB = exp(-2.977 + ln(Ï * D^2 * H))\) |
| \(AGB = Ï * exp(-1.499 + 2.148 * ln(D) + 0.207 * (ln(D))^2 â 0.0281 * (ln(D))^3)\) |
library(BIOMASS)
## For more information on using BIOMASS, visit https://umr-amap.github.io/BIOMASS
## Using temporary cache
## It is recommended to use a permanent cache to avoid to re-download files on each session.
## See function createCache() or BIOMASS.cache option.
woodensity <- getWoodDensity(
genus = NCplot$genus,
species = NCplot$scientificName,
family = NCplot$family
)
## The reference dataset contains 16467 wood density values
## Your taxonomic table contains 93 taxa
head(woodensity)
NCplot$meanWD <- woodensity$meanWD
\(AGB=Ïâexp(â0.667+1.784âln(D)+0.207â(ln(D))2â0.0281â(ln(D))3)\)
NCplot <- NCplot %>%
mutate(dbh_cm = dbh/10,
AGB = meanWD * exp(-0.667 + 1.784*log(dbh_cm) + 0.207*log(dbh_cm)^2 - 0.0281*log(dbh_cm)^3))
(AGB <- sum(NCplot$AGB)/1000) #Mg/ha
## [1] 264.7534
The ratio developed by Mokany et al. (2006). Critical analysis of rootâ:âshoot ratios in terrestrial biomes. Global Change Biology, 12: 84-96 offers specific ratios based on forest type and climate zone and are applicable when the aboveground biomass estimate (shoot) is reported at the stand level
(BGB = 0.275 * AGB)
## [1] 72.80718
Total biomass of 1 hecta
(TB <- AGB + BGB)
## [1] 337.5605
The quantity of carbon can then be estimated by converting biomass to carbon using the IPCC default carbon fraction of 0.47. After that, the carbon stock was multiplied by 44/12 (the carbon atom ratio in the molecular weight of CO2) to get the CO2 equivalent value.
#Carbon stock (ton)
(TC <- TB * 0.47)
## [1] 158.6535
# CO2 (ton)
(CO2 <- TC * 44/22)
## [1] 317.3069
1 carbon credit = 1 ton of CO2 equivalent (tCO2e) emissions.
In 2023, Vietnam successfully sold 10.3 million forest carbon credits for the first time through the World Bank at a price of $5/ton (Global average selling price ~$11.2)
1 hecta forest at NĂși ChĂșa national park ~ 1586 $
(CO2 * 5)
## [1] 1586.535